Non-Makovian decoherence of a two-level system weakly coupled to a bosonic bath 
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Bloch-Redfield equation is a common tool for studying evolution of qubit systems weakly coupled 
to environment. We investigate the accuracy of the Born approximation underlying this equation. 
We find that the high order terms in the perturbative expansion contain accumulating divergences 
that make straightforward Born approximation inappropriate. We develop diagrammatic technique 
to formulate, and solve the improved self-consistent Born approximation. This more accurate treat- 
ment reveals an exponential time dependent prefactor in the non-Markovian contribution dominating 
the qubit long-time relaxation found in Phys. Rev. B 71, 035318 (2005). At the same time, the 
associated dephasing is not affected and is described by the Born-Markov approximation. 



Introduction. A quantum two-level system coupled 
to a bath of harmonic oscillators is a prominent model 
widely employed for describing dissipation in quantum 
physics. The spin-boson model is outlined in several text- 
books [IH3], the advances are covered by the review arti- 
cles (e.g. [MS] and references therein). Immense litera- 
ture is devoted to applications of the model to virtually 
all physics areas ranging from atomic physics and quan- 
tum optics to chemical physics and solid state physics. 
With the advent of quantum information the model re- 
gained increased attention, particularly within the field 
of quantum superconducting circuits [7 . 

For the qubit applications, the most interesting is the 
weak coupling regime when the two-level system slowly 
looses coherence and relaxes to the equilibrium state. 
The Bloch-Redfield equation [51 [U] is a common tool for 
describing this process. It is valid in the lowest, sec- 
ond order approximation with respect to the coupling 
to the bath, and employs the Markov approximation [3]. 
Several schemes beyond these approximations have been 
worked out based on the projection [2,, path integral 
[TJ [5] , diagrammatic [lOj and renormalization group 
techniques. Lifting constraints of the spin-boson model, 
such as a linear coupling to the bath [T^], or equilib- 
rium state of the bath [13l [14], have been discussed in 
literature. However, in the qubit research, the Bloch- 
Redfield equation remains the basic theoretical model, 
whose prediction about an exponential in time qubit evo- 
lution is considered to be qualitatively correct. To what 
extent is this true? The negative answer was obtained 
in [15], where it was found that the qubit long-time re- 
laxation is governed, within the Born approximation, by 
a power-law time dependent non-Markovian term. This 
phenomenon is of a fundamental origin being related to a 
bounded energy spectrum of the bath that generates sin- 
gular branching points of the qubit Green function |16j . 

In this paper we revisit the problem of long-time 
decoherence in the spin-boson model in the weak 
coupling limit. We analyze the perturbation expansion 
and find that the high-order corrections to the Born 
approximation contain accumulating divergences, which 



make straightforward Born approximation inappropriate 
and require summation of the whole perturbation 
series. To identify dangerous perturbative terms and 
perform the summation, we develop a relatively simple 
diagrammatic technique for the Liouville superoperator, 
and formulate an improved, self-consistent Born approx- 
imation. Within this approximation we find that the 
Green function branching points are shifted from the real 
axis yielding an exponential time dependent prefactor 
in the non-Markovian decoherence. For the relaxation, 
the corresponding rate is smaller than the Born rate, 
1/Ti, and thus the non-Markovian term dominates the 
long-time relaxation. For the related dephasing, on 
the contrary, the Born rate, \/2Ti, is larger, and the 
result of the Bloch-Redfield equation remains unchanged. 

Perturbative expansion. The Hamiltonian of the 
spin-boson model has the form 

H = Ha + Hb + Hab, (1) 

where Ha = {ft/2)az is the Hamiltonian of the two- 
level system, Hb = "^k^kb^bkis the Hamiltonian of 
the bosonic bath, and Hat ~ lJ2k''lk {bt + bk) is the 
coupling Hamiltonian, ^ = "fz(Tz + "/^cTx ^ 1 is the cou- 
pling coefficient matrix containing the longitudinal and 
transverse components. We note that in what follows we 
do not use the rotating wave approximation. 

We describe the evolution of the full density matrix, 
p{t), with the Liouville equation in the Laplace form, 

Xp{X)+£p{X) = -tp{0), p{X) = / dte-'^'p{t), (2) 

Jo 

where C is the Liouville (super)operator defined through 
C . . . — [H, . . .], p{0) is the initial condition, which we 
assume in the factorized form, p{0) = Pa{0)pb, with the 
equilibrium density matrix of the bath pt = Ze^^^l"^ . 
The master equation for the reduced density matrix of 
the two-level system, Pa{X) = TTbp{X), is known to have 
the form [2], 

XpaiX) + CaPaW - J{X)pa{X) = -ipa{Q) , (3) 
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FIG. 1. Forth order diagrams for the Green function G: hnes 
represent Go, generally with shifted arguments; dots represent 
vertices R\ dashed lines depict parings. 
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where the self-energy superoperator J{X) represents the 
effect of the bath. Our goal will be to derive a diagram- 
matic representation for the self-energy analogous to the 
one of the many-body diagrammatic theory [17\ 

To this end we introduce the system+bath Green func- 
tion, g{\) through the relation, p{\) — —ig{X)p{0), 
and the Green function, G{X), of the two- level sys- 
tern, pa{X) = -iG(A)pa(0), implying G(A) = (5(A)) = 
TTb[g{X)pb] These Green functions are superoperators 
that satisfy the equations, 

g^9o-9o^ab9, go(A) = (A + £q + (4) 

G = Go + GqJG, GoiX) = {X + Car\ (5) 

Our diagrammatic approach is derived from the observa- 
tion that the n-th order term of the perturbative expan- 
sion of the Green function in Eq. Q , 5 = go ~ go^ab9o + 
. . ., can be exactly presented on the form, 

n / ^ \ 

(A) = ^ J] r/fcA ^+Y1 ^'^fe. ^fe.,..Go(A). (6) 

{k,s}r=l \ i=l ) 

Here we introduced notations, Ru.s ■ ■ ■ — ~ [7&fcs, ■ • ■], 
and 6fc+ = 5^, hk- — bk', the product is ordered, and the 
summation runs over all k and s = ±. To justify this 
equation we consider the action of each term of the ex- 
pansion on some operator of the system, Xa- The action 
of the first term reads, go{X)Xa = Go{X)Xa. The sec- 
ond term follows from the identity, go{X)bksgo{X)Xa — 
Go(A -I- su!k)bksGQ{X)Xa. Continuing this procedure we 
arrive at Eq. and then obtain expansion for G by 
performing the averaging, 

G(A)^5](5(")(A)) (7) 

n 

Averaging in Eqs. ([t]), (|6| is illustrated graphically 
in Fig. [1] The lines represent unperturbed Green func- 
tions. Go, generally with shifted arguments, the dots 
represent superoperatores R, and the dashed lines rep- 
resent the parings of boson operators. All the dia- 
grams split into the two classes: reducible, similar to 
the one in Fig. [l^, and irreducible, Fig. [i}d,c. Since 
the parings select coinciding vertex indices, kr' — kr, 
and Sr' = —Sr, any line of the reducible diagram that 
connects irreducible segments corresponds to Go{X) with 
non-shifted argument. This property also refers to the 
edge lines. Thus any n-th order diagram has the form. 



FIG. 2. Dyson equation and self-consistent self-energy dia- 
grams: bold lines indicate full Green functions G(A), bold 
circle indicates self-energy J{X) containing only crossing di- 
agrams. 



Go(A)(^2™JGo(A)...Go(A)(^2mJGo(A), where (^2™.) 
are irreducible segments, 2mi + • • • + 2mj. ~ n. 

The latter property allows us to formulate an efficient 
diagrammatic theory and present equation for the Green 
function G(A) on the form of the Dyson equation, similar 
to conventional diagrammatics |17j . 

00 

G(A) = Go(A)+Go(A)(5(A))G(A), (5(A)) = ^ (^2™) , 

(8) 

here the summation goes over all irreducible diagrams. 
Comparison of Eqs. ([s]) and ^ establishes explicit form 
for the self-energy, J^{X) — (5(A)). The Born approxi- 
mation corresponds to the first term (^ 7^) in expansion 

Jb(A) - (52) -^r;2(7?fe,_,Go(A-l-sc^fc)i?fc,.). (9) 

k,s 

One may make one step further, and formulate an im- 
proved perturbative expansion by regrouping the series 
terms as shown in Fig. [2] [17] . Here the bold lines (full 
Green functions G) replace the thin lines (free Green 
functions Go) in all the self-energy diagrams, and only 
the diagrams with crossing pairing lines remain. In fact, 
such a rearrangement is necessary because of the presence 
of accumulating divergences in the non-crossing diagrams 
similar to the ones in Fig. [ijj (see below). 



The Dyson equation ([5| with the self-consistent dia- 
grammatic representation for the self-energy in Fig. [2] 
is the main technical result of this paper. The lowest 
order, self-consistent Born approximation is now repre- 
sented with the first diagram in Fig. [2] and is given by 
Eq. ([9]) with the full Green function replacing the free 
Green function, Go — > G. 

For practical calculations, the vector form of the 
density matrix is convenient, then the self-energy be- 
comes the 4x4 matrix. Using the Pauli matrix bases, 
(1, (Tz, f7+, (T_), a± = {ax ± iay)/2, in which the Liou- 
villean of the system acquires the diagonal form, Ca — >■ 



diag(0, 0, il, —Q ), the self-energy matrix takes the form, 
Jb(A)=V/ dujiy{u)P+G{X + suJk) X 

[F+ coth(w/2r) + sP_] . (10) 

Here we introduced the bath spectral density, = 
J2k Vk^i^ ~ ^k)', the matrices P± read, 

p ^ f 7a:P+ \ p = f ^^^^ I'-^P- \ 

+ { 27,p+ j,a, J' - { 27,p„ J ' 

(11) 

where p± = (1/4) [1 =F (^^ + 2(Tzp)], and bar indicates 
conjugated matrix. 



Transverse coupling. We first analyze Eqs. (10), 



(111 for the purely transverse coupling assuming 'jz = 0. 



In this case the P-matrices in Eq. (11) acquire a block- 
antidiagonal form, which induces the block-diagonal form 
of solution to the self-consistent equations ^ and ( [lO| , 
namely, G = diag(G||, Gj_), and Jb = diag(J||, J_i). 
Correspondingly, the master equation for the density ma- 
trix splits into the two independent parts: equation for 
Pz governed by J|| (relaxation) , and equation for p± gov- 
erned by J'± (dephasing). Furthermore, the functional 
dependence between J'b and G acquires peculiar struc- 
ture, 

^|(A) = ^|[Gi(A)], ^i(A) = ^i[G||(A)], (12) 

which allows us to fully investigate the analytical prop- 
erties of the Green functions and thus describe the long- 
time evolution of the density matrix. 

The analytical properties of the Green functions, 
G(A) = (A + £a — Jb(A))^^, are determined by the sin- 
gular points, which consist of the poles and the branch- 
ing points stemming from the self-energy singularities. 
Let us first discuss the positions of the poles; they are 
given by equation, Det(A + £a — Jb(A)) = 0. For the 
longitudinal Green function Gy , the determinant is de- 
generate, and the equation has two roots, A* = 0, and 
A* = J||(A*), where Jy = {J\\)22 is the diagonal element 
of the self-energy. Since Js(A) ^ 7^ is small, the sec- 
ond root is small provided Jy is regular at A = 0. Then 
the approximate position of this pole is Ay = Jy(0), see 
Fig. [3j Similar consideration applies to the transverse 
Green function G_l, whose poles in the main approxima- 
tion are A^ = =F J7 -|- Jj_ (± 17) {J± is the diagonal element 
of JTI). Straightforward evaluation of the pole positions 
using approximation, G ~ Go, in Eq. (10), gives, 



a; = 2mjli.{n) coth ^ = Al = ± r! + ^ (13) 

(in the latter equation, small correction to the real part 
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FIG. 3. Singular points of the longitudinal (left) and trans- 
verse (right) Green functions: poles are indicated with crosses, 
branching points with circles, branch cuts are shown with 
dashed lines. 



The branching points appear in the self-energy as the 
result of integration over uj in Eq. ( 10 ) of the singular 



Green function: Since the integral is defined on the semi- 
axis, w > 0, the singularity builds up when the Green 
function pole approaches the integration edge a; = 0. 
Consequently, the branching points, Ay, of J'\\ coincide, 
by virtue of Eq. (12), with the poles of the transverse 
Green functions. Ay = A^, and vice versa for the trans- 
verse part, X'^ = Ay as shown on Fig. jSj This veri- 
fies the above made assumption regarding analytical be- 
havior of the self-energy near the poles of corresponding 
Green function. 

Important conclusion of our analysis is that both the 
poles and branching points lie in the complex A-plane. 
Therefore the time evolution of the density matrix is al- 
ways exponential. This conclusion is the main result of 
the paper. 

Within the straightforward Born approximation of 
Eq. (|9|, the branching points lie on the real axis be- 
ing defined by the poles of non-perturbed Green func- 
tions, Gq. This leads to the non-exponential time re- 
laxation found in |15j . However, the higher order dia- 
grams also contain the singularities that accumulate with 
the growing diagram order. The most dangerous dia- 
grams are of the type of Fig. [TJ), having the singularities, 
7k™/ ^ A*)™. This shortcoming of the straightforward 
perturbation theory is eliminated by the improved per- 
turbative expansion. 

A more accurate treatment within the self-consistent 
Born approximation reveals the shifts of the self-energy 
branching points into the complex plane. In the case of 
Gy, they situate more close to the real axis compared 
to the pole, left panel in Fig. |3j and thus define the the 
long-time relaxation. Explicit solution for Pzit) has the 
form. 
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The contribution of the poles reads, 



1 n ^ 

- tann — - e 

2 2T^ 



-rt 



1). (15) 



This coincides with the result of the Markovian approx- 
imation of the Bloch-Redfield theory 9 , T = l/Ti. The 
result of integration along the branch cuts depends on 
the specific form of the bath spectral density. Here we 
assume the ohmic bath for certainty, i'{uj) = Kijje~^/^'' ^ 
with a large cutoff frequency, ^ S> T. The asymp- 
totical contribution of the branch cuts at i7t ^ 1 reads, 



PI 



where 



F{T,t) 



2T 



fi2Tt), 



sinh^ {-KTt) 
f{x)^x^ / e~^^y^ cot ?/. 



(16) 



(17) 



The function F{T,t) approaches the asymptotic values, 
F = 1 in the limit Tt < 1, and F = 2T/ojc at Tt > 1. 



Although the contribution from branch cuts in Eq. ( 16 ) 
decays slower than the Markovian part, Eq. (151, it 



contains a small rapidly oscillating factor, and because 
of this, the non-Markovian part becomes dominant at 
rather large times. 



t > 4Ti In — > Ti 



(18) 



The oscillation, oc cos fit, formally results from the po- 
sitions of the branching points, defined by the poles of 
Gj_, away from the imaginary axis. It can therefore be 
interpreted as the effect of coupling to the transverse evo- 
lution (cf. Eq. ([l9| below). 

The situation is different for the dephasing. In this 
case the branching point of G±, right panel in Fig. |3j 
has larger imaginary part than the poles, and thus non- 
Markovian contribution decays more rapidly. Therefore 
the long-time dephasing in the case of transverse coupling 
is given by the Markov approximation, 



pY-^(t)=p+(0)e-^"*-(r/2)*, 



'+ 



^ ' (m)2 



(p+(0) + p_(0)). (19) 



General coupling. The presence of both the trans- 
verse and longitudinal couplings, Iz 7^ Oj in Eq- (H) 



leads to the two complications: First, the block-diagonal 
form of the Green functions is lost due to the fact that 
the P-matrices acquire diagonal parts, and the relax- 
ation and dephasing become mixed. Second, the prop- 



erty, Eq. (12), is lost, which implies more complex struc- 



ture of the singularities consisting of the coinciding poles 
and branching points. Attentive analysis, however, shows 



that the mixing of the relaxation and dephasing is of the 
higher order, ~ 7**, and within the Born approximation, 
this effect can be neglected. Furthermore, the property, 
Eq. (12 1, is only violated for the transverse (dephasing) 



part of the self-energy, 



-J'AGi 



(20) 



while for the longitudinal part it persists. Thus, up to the 
fourth-order corrections, the relaxation is not affected by 
the presence of z-coupling. 

For the dephasing, the major problem is the contri- 
bution of the singular points at A = ifi -|- ir/2, af- 
fected by the second term in Eq. (20). To evaluate this 
contribution we note that the first term in Eq. ( 20 ) , 



J7_l[G||] oc 72;, is related to the transverse coupling, and it 
is analytical at this point and can be evaluated as before, 
Jj_[G||] w Jl(±f7) = ir/2. Thus the role of this term is 
to shift the frequency, — >■ cr^ri — ir/2, in the equation 
for the transverse Green function, 



Gi(A) = (A + a,^l - iT/2 - J'_^{\)Y 



(21) 



This equation, however, apart from the frequency shift, 
describes the effect of pure z— couphng (since J'^ oc jz), 
and the solution to this case is well known in the lit- 
erature (e.g. [7] and references therein). Therefore the 
dephasing in the present case is given by this solution 
with the shifted frequency. 



p+(t)^-»^^*-(r/2)ti^(t) p^(o) 
F(i)=exp 



i/(w) 
acu — ^ sm 

^ 2 



2'"*coth^ 



(22) 



Conclusion. We investigated long-time decoherence 
of a two-level system weakly coupled to a bosonic bath. 
We found that the conventional Born approximation does 
not correctly describe the evolution of the system, and 
formulated an improved, self-consistent Born approxima- 
tion based on a developed diagrammatic technique. We 
found an exponential time dependent prefactor in slow 
varying non-Markovian tail that dominates the long-time 
relaxation, while the dephasing associated with relax- 
ation is well described within the Markov approximation. 
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